
function [model_moments_bymkt, pareto_weight_evol, U_ud_bymkt, U_ud_tot_bymkt] =  moments_ud_DivorceSettlement(n_g, T, lt, dt, H_max, disc, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, home_prod, match_qual, m_th, n_th, gamma, L_max, params, lambda_ud, data_f, data_m);

parest(1,:)=params(1:9); 
parest(2,:)=params(10:18); 
sigma_f=[params(19) params(19) params(19) params(20) params(20) params(20) params(21) params(21) params(21)]; 
sigma_m=[params(22) params(23) params(24) params(22) params(23) params(24) params(22) params(23) params(24)]; 
match_qual(:,2)=params(25:33)'; 

%- Wage rate
a0m_1=params(34);
a0m_2=params(35);
a0m_3=params(36);

a0f_1=params(37);
a0f_2=params(38);
a0f_3=params(39);

%- Returns to experience
a1m_1=params(40);
a1m_2=params(41);
a1m_3=params(42);

a1f_1=params(43) ;
a1f_2=params(44) ;
a1f_3=params(45) ;

%- Returns to experience^2 and beyond
a2m_1=params(46);
a3m_1=0;

a2m_2=params(47);
a3m_2=0;

a2m_3=params(48);
a3m_3=0;

a2f_1=params(49);
a3f_1=0;

a2f_2=params(50);
a3f_2=0;

a2f_3=params(51);
a3f_3=0;

%- Returns to sahw
bm_1=params(52);
b2m_1=0;
bm_2=params(53);
b2m_2=0;
bm_3=params(54);

b2m_3=0;
bf_1=0;
b2f_1=0;
bf_2=0;
b2f_2=0;
bf_3=0;
b2f_3=0;

det_earn_m_1=[a0m_1, a1m_1, a2m_1, a3m_1, bm_1, b2m_1];
det_earn_m_2=[a0m_2, a1m_2, a2m_2, a3m_2, bm_2, b2m_2];
det_earn_m_3=[a0m_3, a1m_3, a2m_3, a3m_3, bm_3, b2m_3];

det_earn_m=[det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3];

det_earn_f_1=[a0f_1, a1f_1, a2f_1, a3f_1, bf_1, b2f_1];
det_earn_f_2=[a0f_2, a1f_2, a2f_2, a3f_2, bf_2, b2f_2];
det_earn_f_3=[a0f_3, a1f_3, a2f_3, a3f_3, bf_3, b2f_3];

det_earn_f=[det_earn_f_1;det_earn_f_1;det_earn_f_1;det_earn_f_2;det_earn_f_2;det_earn_f_2;det_earn_f_3;det_earn_f_3;det_earn_f_3];

parest(3,:)=lambda_ud; 

[dec, lambda_grid, pareto_weight_index, U_f, U_m, U_f0, U_m0, U_ftot, U_mtot] = solution_backwards_ud_DivorceSettlement(n_g, T, lt, dt, H_max, disc, det_earn_f, det_earn_m, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, parest(1,:), home_prod, parest(2,:), match_qual, m_th, n_th, gamma, L_max, parest(3,:), sigma_f, sigma_m);

I = 10000;

hw = ones(n_g,I,T+1); 
lc = zeros(n_g,I,T+1); 
lcd=zeros(n_g,I,T-lt+1); 
hwd=zeros(n_g,I,T-lt+1); 

n_p=4;

married_moments=zeros(n_g,n_p); 
singles_moments=zeros(6,1); 
for g=1:n_g
married_moments(g,1)=(exp(U_f(g))/U_ftot(g))*data_f(g); 
married_moments(g,2)=(exp(U_m(g))/U_mtot(g))*data_m(g); 
end

singles_moments(1,1)=(exp(U_f0(1))/U_ftot(1)); 
singles_moments(2,1)=(exp(U_f0(4))/U_ftot(4));
singles_moments(3,1)=(exp(U_f0(7))/U_ftot(7));

singles_moments(4,1)=(exp(U_m0(1))/U_mtot(1)); 
singles_moments(5,1)=(exp(U_m0(2))/U_mtot(2));
singles_moments(6,1)=(exp(U_m0(3))/U_mtot(3));

rp=zeros(n_g,I,T+1);
for g=1:n_g
   [~, rp(g,:,:)]=min(abs(lambda_grid-parest(3,g)));  
end

pareto_weight=NaN(n_g,I,T-lt+1);
pareto_weight_evol=ones(n_g,T);

%......................... Simulate ........................% 
seed=1812;
rng(seed); 

wft_rnd_all=rand(3,I,T); 
wft_rnd_all_rep=[wft_rnd_all(1,:,:) ; wft_rnd_all(1,:,:) ;wft_rnd_all(1,:,:) ;wft_rnd_all(2,:,:) ; wft_rnd_all(2,:,:) ;wft_rnd_all(2,:,:); wft_rnd_all(3,:,:) ; wft_rnd_all(3,:,:) ;wft_rnd_all(3,:,:) ];
wmt_rnd_all=rand(3,I,T); 
wmt_rnd_all_rep=[wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ];
tht_rnd_all=rand(9,I,T); 
wft_rnd=zeros(I,T); 
wmt_rnd=zeros(I,T); 
tht_rnd=zeros(I,T); 
wft_gridpt=zeros(I,T);
wmt_gridpt=zeros(I,T);
th_gridpt=zeros(I,T);

for g=1:n_g 
%-Female wages
[~,~,cdf_fu,cdf_fc]=wageproc_fn(stoch_earn_f(g,:), m_eta_f, n_eta_f);
%-Male wages
[~,~,cdf_mu,cdf_mc]=wageproc_fn(stoch_earn_m(g,:), m_eta_m, n_eta_m);
%-Theta process
[~,~,cdf_thu,cdf_thc]=thetaproc(match_qual(g,:), m_th, n_th);

wft_rnd(:,:)=wft_rnd_all_rep(g,:,:); 
wmt_rnd(:,:)=wmt_rnd_all_rep(g,:,:); 
tht_rnd(:,:)=tht_rnd_all(g,:,:);

%- Female wage nodes
for i=1:I   
       for t=lt:T 
           if t==lt; 
        for f=1:n_eta_f 
            if wft_rnd(i,t)<=(cdf_fu(1,f));
               wft_gridpt(i,t)=f;
               pastnode=f;
               break
            end
        end
        
           else 
               
        for f=1:n_eta_f 
            if wft_rnd(i,t)<= (cdf_fc(pastnode,f));
            wft_gridpt(i,t)=f;
            newnode=f;
            break
            end
        end
        
        pastnode=newnode;
                   
           end
           
       end
end

%- Male wage nodes
for i=1:I   
       for t=lt:T 
            
           if t==lt;
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<=(cdf_mu(1,m));
               wmt_gridpt(i,t)=m;
               pastnode=m;
               break
            end
        end
        
           else
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<= (cdf_mc(pastnode,m));
            wmt_gridpt(i,t)=m;
            newnode=m;
            break
            end
        end
        pastnode=newnode; 
        
           end
           
       end
end

%- Theta nodes

for i=1:I   
       for t=lt:T 
           if t==lt;
        for th=1:n_th 
            if tht_rnd(i,t)<=(cdf_thu(1,th));
               th_gridpt(i,t)=th;
               pastnode=th;
               break
            end
        end
        
           else
               
        for th=1:n_th  
            if tht_rnd(i,t)<= (cdf_thc(pastnode,th));
            th_gridpt(i,t)=th;
            newnode=th;
            break
            end
        end
        
        pastnode=newnode; 
        
           end
           
       end
end


%......... Compute choices ..................%

hw(g,:,lt+1:T+1) = 0; 

for i=1:I
   
    for t=lt 
        lc(g,i,t) = dec(g, rp(g,i,t) ,hw(g,i,t), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t);
        if lc(g,i,t) == 2
            hw(g,i,t+1) = hw(g,i,t) + 1; 
        elseif lc(g,i,t) == 1
            hw(g,i,t+1) = hw(g,i,t);
        end
        
        rp(g,i,t+1)= pareto_weight_index(g,rp(g,i,t),hw(g,i,t),wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t),t); 
        pareto_weight(g,i,t)=lambda_grid(rp(g,i,t+1));
         end
      
    for t=lt+1:T 
        if lc(g,i,t-1)==3 
            lc(g,i,t)=-1;            
        elseif lc(g,i,t-1)==-1 
            lc(g,i,t)=-1;
            
        else 
        lc(g,i,t) = dec(g, rp(g,i,t), hw(g,i,t), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t); 
        end
        if lc(g,i,t) == 2
            hw(g,i,t+1) = hw(g,i,t) + 1; 
            rp(g,i,t+1)= pareto_weight_index(g,rp(g,i,t),hw(g,i,t),wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t),t);             
            pareto_weight(g,i,t)=lambda_grid(rp(g,i,t+1));                                                                                                       
        elseif lc(g,i,t) == 1
            hw(g,i,t+1) = hw(g,i,t);
            rp(g,i,t+1)= pareto_weight_index(g,rp(g,i,t),hw(g,i,t),wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t),t);
            pareto_weight(g,i,t)=lambda_grid(rp(g,i,t+1));  
        else
            hw(g,i,t+1) = 0; 
            rp(g,i,t+1)=NaN; 
        end
        
        
    end
end

lcd(g,:,:)=lc(g,:,lt:T); 
hwd(g,:,:)=hw(g,:,lt:T);

%......... Moments ..................%
fracd=lcd==3;
frachw=lcd==2;
divorced=(lcd==3 | lcd==-1);
married_moments(g,3)=sum(sum(frachw(g,:,1:T)))/(I*T-sum(sum(divorced(g,:,1:T)))); 
married_moments(g,4)=sum(sum(fracd(g,:,1:T)))/I; 

for t=lt:T
pareto_weight_evol(g,t)=nanmean(pareto_weight(g,:,t));
end

end 

model_moments_bymkt = [married_moments(:,1); married_moments(:,2); singles_moments; married_moments(:,3); married_moments(:,4)];
U_ud_bymkt=[U_f U_f0 U_m U_m0]; 
U_ud_tot_bymkt=[U_ftot U_mtot]; 

end


